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Abstract.  Molecular  dynamics  simulations  of  vapor-liquid  equilibrium  states  and  those  of  evaporation  from  liquid  phase  into 
a  virtual  vacuum  are  performed  for  water.  In  spite  of  the  formation  of  molecular  clusters  in  the  vapor  phase  and  the  presence 
of  the  preferential  orientation  of  molecules  at  the  interface  due  to  uneven  sharing  of  the  bonding  electron  pair,  essentially  the 
same  results  as  in  our  previous  study  for  argon  are  obtained.  That  is,  when  the  bulk  liquid  temperature  is  relatively  low,  the 
distribution  function  of  evaporation  can  be  expressed  as  the  product  of  the  equilibrium  distribution  of  saturated  vapor  at  the 
temperature  in  the  bulk  liquid  phase  and  a  well-defined  evaporation  coefficient,  which  is  determined  as  a  decreasing  function 
of  the  liquid  temperature,  and  is  found  to  approach  unity  with  the  decrease  of  the  temperature. 


INTRODUCTION 


The  kinetic  boundary  condition  at  an  interface  between  a  vapor  and  its  condensed  phase 


/out  =  af  +  (1  -  a)f  (4  >  0) , 


(1) 


has  widely  been  used  so  far,  where  /out  is  the  distribution  function  of  outgoing  molecules  from  the  interface,  a  is  a 
parameter  between  zero  and  unity,  sometimes  called  the  condensation  coefficient,  /e  is  the  equilibrium  distribution 
of  saturated  vapor  at  the  temperature  of  the  condensed  phase,  f  is,  usually,  the  distribution  function  of  the  diffuse 
reflection  at  the  temperature,  and  qz  is  the  velocity  component  normal  to  the  interface  (for  simplicity,  we  consider  a 
planar  interface  at  rest,  facing  the  positive  z  direction).  We  call  a/c  the  evaporation  part  and  (1  —  a)/r  the  reflection 
part. 

Many  important  phenomena  associated  with  evaporation  and  condensation  at  the  interface  have  been  clarified  by 
the  use  of  this  type  of  kinetic  boundary  condition  (see  [1,  2,  3]  and  references  therein).  However,  its  physical  validity 
still  remains  to  be  established,  and  the  related  studies  on  the  basis  of  the  molecular  dynamics  (MD)  method  have  just 
been  started  [4,  5,  6,  7].  In  a  previous  paper  [7],  to  study  the  kinetic  boundary  condition  at  an  interface  between  argon 
vapor  and  its  condensed  phase,  we  have  carried  out  MD  simulations  of  vapor-liquid  (or  solid)  equilibrium  states  (see 
Fig.  1)  and  those  of  evaporation  from  the  liquid  (or  solid)  phase  into  a  virtual  vacuum.  As  the  result,  we  have  shown 
that,  in  the  case  that  the  temperature  of  the  bulk  condensed  phase  7)  is  near  the  triple  point  temperature  of  argon, 
the  distribution  function  for  molecules  evaporating  into  vacuum  is  given  by  aepvf* ,  where  pv  is  the  saturated  vapor 
density,  pvf*  is  the  half-Maxwellian  of  saturated  vapor 


Pvf*  = 


{InRTffi!2 


exp  | 


S2 


2  RTp 


)  (4>0), 


(2) 


(R  is  the  gas  constant  and  qx  and  are  the  velocity  components  tangential  to  the  interface),  and  ae  is  the  evaporation 
coefficient  defined  by  the  ratio  of  an  averaged  mass  flux  evaporating  into  vacuum  (Jev ap)  to  averaged  outgoing  mass 
flux  in  the  equilibrium  state  (/out}eq> 


ae 


(3) 


where  the  brackets  denote  the  ensemble  average  (see  Fig.  2). 
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FIGURE  1.  Water  vapor-water  equilibrium  in  a 
simulation  cell. 
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FIGURE  2.  Averaged  mass  fluxes  at  the  vapor-liquid  interface. 


In  this  paper,  we  shall  extend  the  previous  study  to  the  case  of  water.  The  simulations  are  carried  out  in  a  similar 
way  to  the  previous  one.  However,  a  more  careful  and  detailed  analysis  is  required  because  a  water  molecule  has  the 
additional  degrees  of  freedom  for  internal  motions  and  the  polarity  due  to  uneven  sharing  of  the  bonding  electron 
pair.  In  particular,  the  latter  induces  the  formation  of  molecular  clusters  in  vapor  phase  and  leads  to  the  presence 
of  the  preferential  orientation  of  molecules  at  the  interface.  Nevertheless,  as  shown  in  the  present  study,  we  arrive 
at  essentially  the  same  conclusions  as  the  previous  study.  That  is,  the  distribution  function  of  evaporation  can  be 
expressed  as 

0CePvf*8*  (&>0),  (4) 

in  relatively  low  temperature  cases,  where  pvf*  is  the  half-Maxwellian  for  translational  motions  of  center  of  mass  of  a 
polyatomic  molecule,  which  is  equal  to  that  given  in  the  right-hand  side  of  Eq.  (2)  for  a  monatomic  molecule,  and  g* 
is  the  equilibrium  distribution  associated  with  internal  motions  of  polyatomic  molecule  of  n  degree  of  freedom. 


£"/ 2~1  /  E  \ 

E(«/2)  (kTe)n/2  6XP  \  kTj 


(0  <E  <°°), 


(5) 


(k  is  the  Boltzmann  constant,  E  is  the  energy  of  internal  motion  of  one  molecule,  and  T  is  the  gamma  function).  The 
symbol  ~  signifies  a  normalized  distribution  function  and  the  superscript  *  represents  that  the  distribution  of  molecules 
is  in  an  equilibrium  state. 


METHOD  OF  ANALYSIS 

As  in  the  previous  study  [7],  we  first  define  a  distribution  function  of  spontaneous  evaporation,  /ev ap,  as  a  distribution 
function  of  molecules  evaporating  from  the  interface  and  independent  of  the  condition  of  incident  vapor  molecules. 
This  means  that  /eva p  is  unchanged  whether  the  condensed  phase  is  contact  with  the  vapor  or  exposed  to  vacuum  and  it 
should  be  determined  by  the  temperature  in  the  bulk  condensed  phase  only.  An  arbitrary  /out  given  by  MD  simulation 
can  then  be  split  into  two  parts, 

/out  =  /ev ap  +  /ref  (4  >  0) ,  (6) 

where  /ref  =  /out  —  ,/evap  ■  Putting  /out  in  the  form  of  Eq.  (6)  enables  us  to  verify  the  evaporation  part  of  Eq.  (1)  in 
the  method  explained  below.  The  splitting  may  also  be  regarded  as  the  extraction  of  an  inherent  property  of  the  bulk 
condensed  phase  from  /out.  We  shall  remark  that  we  don’t  intend  classifying  every  individual  outgoing  molecule  as 
either  evaporated  or  reflected  one  in  a  vapor-liquid  two-phase  system;  the  complete  classification  may  be  impossible 
[6,7]. 

Suppose  that  there  are  no  incoming  and  reflected  molecules  at  the  interface,  and  then  we  have  from  Eq.  (6) 

/out  =  /evap  (4>0)-  (7) 

Such  an  extreme  situation  can  be  realized  in  an  MD  simulation  by  eliminating  vapor  molecules  at  a  distance  near  the 
interface.  We  here  call  this  the  evaporation  into  virtual  vacuum,  or  simply,  evaporation  into  vacuum.  The  functional 


492 


FIGURE  3.  Saturated  vapor  density. 


FIGURE  4.  Boundary  condition  and  temperature  control. 


form  of  /evap  can  therefore  be  determined  in  the  MD  simulation  of  evaporation  into  vacuum  by  counting  the  number 
of  evaporating  molecules  and  sampling  its  velocity. 


NUMERICAL  METHOD 


Intermolecular  potential 


As  an  intermolecular  potential,  we  adopt  TIP3P  model  for  water  [8],  since  it  has  been  widely  used  and  the 
computation  time  is  not  so  large;  we  have  to  run  many  simulations  to  obtain  a  large  number  of  samples  for  the  accurate 
construction  of  the  distribution  function.  The  potential  is  a  3-site  rigid  model  consist  of  Lennard-Jones  12-6  potential 
and  Coulomb  potential.  The  additional  degree  of  freedom  is  n  =  3  in  Eq.  (5)  and  E  is  the  internal  energy  associated 
with  the  rotation  of  molecule.  The  intermolecular  potential  between  a  site  m  of  a  molecule  and  a  site  n  of  a  different 
molecule  can  be  written  as 


tymn  —  4£n 


f 

(  Gmn 

1 

1  1  ^mn 

2 

\  ? mn  J 

\  r mn  J 

4neo 

rmn  r3ut 

r cut 

(8) 


where  rmn  is  the  distance  between  the  site  m  and  the  site  n,  z„,e o  and  zne o  are  partial  charges  for  the  sites,  eo  is  the 
elementary  charge,  £o  is  the  dielectric  constant  of  vacuum,  and  emn  and  omn  are  the  usual  Lennard-Jones  parameters. 
As  indicated  in  Eq.  (8),  the  electrostatic  term  is  shifted  and  scaled  smoothly  to  zero  at  rcut  (0.9  nm),  and  lattice 
summation  techniques  are  not  used.  The  applicability  of  the  shifted  and  scaled  Coulomb  potential  has  been  confirmed 
in  [9],  and  it  is  known  that  the  use  of  the  Ewald  summation  to  the  vapor-liquid  two-phase  system  leads  to  some 
unnatural  behaviors  [10].  The  Lennard-Jones  part  also  is  truncated  at  the  same  cutoff  distance  rcut. 


Vapor-liquid  equilibrium 

In  the  equilibrium  simulations,  the  periodic  boundary  conditions  are  imposed  for  all  three  directions  of  simulation 
cell.  In  lower  temperature  cases  of  T(  <  350  K,  we  use  the  cell  of  50  x  50  x  200  A3,  in  which  N  =  2000  molecules 
are  contained,  and  at  higher  temperatures  7>  >  360  K,  we  treat  4000  molecules  in  the  longer  cell  of  50  x  50  x 
400  A.  After  equilibrating  the  system,  we  continue  the  simulation  for  5  ns  and  accumulate  the  configurations  of 
all  molecules  in  the  cell  every  1  ps.  Since  the  system  has  two  interfaces  as  shown  in  Fig.  1,  the  number  of  samples 
of  configurations  is  Ns  =  10000.  The  ensemble  averages  for  various  macroscopic  quantities  can  be  evaluated  from  Ns 
sampled  configurations. 

At  lower  temperatures,  the  saturated  vapor  densities  calculated  with  TIP3P  model  almost  agree  with  experimental 
values.  However,  the  calculated  value  of  saturated  vapor  density  becomes  large  compared  with  experimental  one  as  the 
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FIGURE  5.  Temporal  evolution  of  averaged  mass  flux,  (a):  T*  =  0.60  (Tt  =  340  K);  (b):  T*  =  0.72  (Te  =  420  K). 


FIGURE  6.  Spacial  distribution  of  averaged  mass  flux. 


temperature  increases.  Therefore,  the  present  results  are  not  directly  compared  with  experimental  ones.  To  compensate 
the  disagreements  of  saturated  vapor  density,  we  arrange  the  results  in  terms  of  a  reduced  temperature  T*  defined  by 

y,*  _  -^triple  T  A(7^  ^triple) 

Ter 

where  A  is  a  fitting  parameter  and  7’tripie(=  273.2  K)  and  rcr(=  647.1  K)  are  experimentally  evaluated  triple  point  and 
critical  point  temperatures,  respectively  [1 1].  As  shown  in  Fig.  3,  A  =  4/3  gives  good  agreement  for  water. 


Evaporation  into  virtual  vacuum 


In  the  vacuum  simulations,  we  introduce  an  open  boundary  at  a  distance  from  the  interface  and  eliminate  molecules 
there  (Fig.  4),  while  the  periodic  boundary  conditions  are  applied  in  the  x  and  y  directions.  Molecules  evaporate  into 
the  virtual  vacuum  through  the  open  boundary  and  the  interface  recedes  with  time  as  a  result  of  the  evaporation.  The 
steady  evaporation  state  is  realized  on  the  moving  coordinate 


*  Z  (Zm  )  /v 

Z  —  5  )  Vs  —  1 

O  pe 


(10) 


where  <5  is  the  10-90  thickness  of  the  transition  layer  whose  center  is  located  at  Z,„  on  the  original  z  coordinate,  t  is  the 
time  from  the  beginning  of  the  vacuum  simulation,  Js  is  a  molecular  flux  evaporating  into  vacuum,  and  vs  is  the  speed 
of  the  moving  coordinate.  The  averaged  values  for  macroscopic  quantities  and  the  distribution  functions  are  evaluated 
on  the  moving  coordinate  z*. 

The  control  of  the  temperature  in  the  bulk  liquid  phase  is  essential  for  the  realization  of  the  steady  evaporation  state. 
Using  the  velocity  scaling  method,  we  control  the  temperature  of  the  liquid  phase  in  the  region  z*  <  —L*  as  shown 
schematically  in  Fig.  4.  As  the  result,  the  averaged  temperature  of  the  bulk  liquid  phase  is  kept  almost  uniform  and 
constant  at  a  specified  7}-.  See  Ref.  [7]  for  the  detail  of  the  simulation  technique. 
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FIGURE  7.  Distribution  function  of  translational  velocity  of  the  center  of  mass  of  molecule. 


FIGURE  8.  Distribution  function  of  angular  velocity  around  the  principal  axis  of  water  molecule.  The  principal  axes  are  signified 
by  x,  y,  and  z. 


RESULTS  AND  DISCUSSION 
Evaporation  mass  flux 

We  shall  start  with  the  verification  of  steadiness  of  the  evaporation  into  vacuum,  because  the  realization  of  the  steady 
evaporation  state  implies  the  existence  of  the  spontaneous-evaporation  mass  flux  (yevap)  determined  by  the  bulk  liquid 
temperature  7).  only.  After  an  initial  transient  effect  is  died  out,  almost  steady  evaporation  state  is  established  as  shown 
in  Fig.  5.  The  spatial  uniformity  of  the  net  mass  flux  is  clearly  shown  in  Fig.  6. 


Distribution  functions 

Now,  we  shall  evaluate  the  distribution  function  of  the  translational  velocity  of  molecules  evaporating  into  vacuum, 
/trans,  and  that  of  internal  energy  associated  with  rotational  motion,  grot-  In  Fig.  7,  the  velocity  distributions  of 
evaporating  molecules  evaluated  at  z*  =  L*  are  plotted  for  some  temperatures.  The  abscissa  £)  =  | j/y/lRT^  is  the  j 
component  of  the  normalized  molecular  velocity.  Figure  7  shows  that  the  translational  velocity  distributions  of  £Y  and 
Cv  denoted  by  triangles  and  squares  almost  agree  with  a  one-dimensional  normalized  Maxwellian  (l/v/7r)exp(— £?) 
denoted  by  a  solid  curve,  although  those  for  relatively  high  7)  cases  slightly  shift  to  lower  temperature  distributions. 
For  relatively  low  7)  cases  shown  in  Figs.  7(a),  the  velocity  distribution  of  denoted  by  closed  circles  becomes  nearly 
a  one-dimensional  normalized  half-Maxwellian  (2/\Jn)  exp(— £?)  (£z  >  0)  denoted  by  a  dashed  curve. 

In  Fig.  8,  we  plot  the  angular  velocity  distributions  of  evaporating  molecules  at  z*  =  L*  The  abscissa  Vj  = 
0)j  \/ lj / (2kTf)  in  the  figure  denotes  the  normalized  angular  velocity  component  around  the  principal  axis.  As  can  be 
seen  from  Fig.  8,  the  angular  velocity  distributions  are  isotropic  and  are  nearly  the  Maxwellians  for  all  temperatures, 
although  the  vapor  is  in  an  extreme  nonequilibrium  condition.  As  in  the  case  of  ftrws,  grot  also  slightly  shift  to  lower 
temperature  distributions  in  higher  temperature  cases. 

Consequently,  the  distribution  function  of  molecules  evaporating  into  vacuum  /evap  may  be  written  as  the  product 

of  Pc,  /trans>  alld  <?rot’ 

/evap  =  Pc  /trans  .grot  =  ~Pcf  /  ■  (11) 
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FIGURE  9.  Evaporation  and  outgoing  mass  fluxes. 


FIGURE  10.  Evaporation  coefficient. 


where  pc  is  the  density  of  the  vapor  evaporating  into  vacuum.  Using  the  fact  that  /trans  is  the  normalized  half- 
Maxwellian,  it  is  easy  to  show  that  pc  =  \(Xepv  [7],  and  as  a  result,  we  obtain 

/evap  =  aePvf*f  =  aef.  (12) 

The  evaporation  coefficient  as  a  function  of  T*  is  shown  in  Fig.  10,  where  the  previous  results  for  argon  [7]  and  for 
methanol  [12]  are  also  plotted. 


CONCLUSIONS 

We  have  carried  out  the  MD  simulations  of  vapor-liquid  equilibrium  and  those  of  steady  evaporation  into  the  virtual 
vacuum  for  water.  The  distribution  function  of  molecules  evaporating  into  vacuum  has  been  accurately  obtained.  We 
have  demonstrated  that  in  relatively  low  temperature  case  the  distribution  function  is  the  product  of  the  evaporation 
coefficient,  the  half-Maxwellian  of  translational  molecular  velocity,  and  the  equilibrium  distribution  of  rotational 
energy.  The  evaporation  coefficients  of  water  is  also  determined  as  a  decreasing  function  of  the  bulk  liquid  temperature, 
and  their  values  are  found  to  become  close  to  unity  with  decrease  in  the  temperature. 
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